%% initialize
clear all;
close all;
clc

cd  E:\ReplicateBuild\02_code\03_matlab


% diary('mixed_logit')
% diary off
rng(328120)

% load principal prefs and clean them up

  load '../../01_data/05_temp/mixed_logit_principal_totVA_extraX' 
    
    IJX_princ=[ ones(length(IJX_princ),1) IJX_princ];   % add a vector of ones  
    IJ_princ=IJX_princ * thetahat(1:end-1,1);     % don't include RC
    clearvars -except IJ_princ IJN* IJVA* exper_potential

  load('../../01_data/05_temp/mixed_logit_rc6_totVA_DISAD_Ndis')   
  
   savefile= '../../01_data/05_temp/model/thetahat';
  save(savefile,'thetahat');
 
    
    % select the best Xdraws and clean them up
    [Xdraws_best] = func_likely_draws_current(IJitidx, IJchoiceset,IJcurrentFirst, Xdraws,Xmat,thetahat,y_logistic);
    
    % redo Xdraws_best to be as long as everything else
    [Ci, IAi, ICi]=unique(IJitidx(IJchoiceset==1));
    Xdraws_unique=Xdraws_best(IAi,:);
    Xdraws_best2=Xdraws_unique(IJitidx,:);
    
    % find predicted values in each match (teacher preferences)
    
    IJipref=[Xmat Xdraws_best2]*thetahat;
    
    sigma=std(IJipref);
    mu=mean(IJipref);
    
    IJipref_rand=normrnd(mu,sigma,length(IJipref),1);
    
    % compare variance of preferences, and then simulate them randomly
    
    % generate IJpref that is -999999999 if choiceset=0
    
    IJipref_choice=IJchoiceset.*IJipref + (1-IJchoiceset)*(-9999999999);
    
    IJipref_rand_choice=IJchoiceset.*IJipref_rand + (1-IJchoiceset)*(-9999999999);

    
    % omit the RCs
    IJipref_norc= Xmat(:,1:36)*thetahat(1:36,1);
    
   IJipref_norc_choice=IJchoiceset.*IJipref_norc + (1-IJchoiceset)*(-9999999999);
     
    % build school preferences (by changing this block, can make it principal preferences, or....)
      
    IJoutput=IJclassSize.*IJVAma;
    IJoutput_dis=IJclassSize.*IJfracDisadv.*IJVAm2;
    IJoutput_adv=IJclassSize.*(1-IJfracDisadv).*IJVAm1;
    IJN_dis=IJclassSize.*IJfracDisadv;
    
  % generate ouput that is -999999999 if choiceset=0
    
    IJoutput_choice=IJchoiceset.*IJoutput + (1-IJchoiceset)*(-9999999999);   
% generate IJ_princ that is -99999 if choiceset=0

    IJ_princ_choice=IJchoiceset.*IJ_princ + (1-IJchoiceset)*(-9999999999);  
    IJ_choice=(1-IJchoiceset)*(-9999999999);
    IJ=zeros(length(IJchoiceset),1);
    
        % generate IJfracDisadv that is -99999999 if choicset=0
    
    IJfracDisadv_choice=IJchoiceset.*IJfracDisadv + (1-IJchoiceset)*(-9999999999);
 
       % generate N disadvantage
       
           IJNdis= IJclassSize.*IJfracDisadv;
    IJNdis_choice=IJchoiceset.*IJNdis + (1-IJchoiceset)*(-9999999999);
    
    % generate principal for non-title I, and VA for title I
    
    IJ_princ_ntitleI_va_titleI=IJ_princ.*(IJfracDisadv<0.4) + IJoutput.*(IJfracDisadv>0.4)+IJoutput.*(IJfracDisadv==0.4);
    IJ_princ_titleI_va_ntitleI=IJ_princ.*(IJfracDisadv>0.4) + IJoutput.*(IJfracDisadv<0.4)+IJ_princ.*(IJfracDisadv==0.4);

    % generate 4 variables from the x matrix
    
    IJblack=Xmat(:,30);
    IJexp23=Xmat(:,33);
    IJexp46=Xmat(:,34);
    IJexp7p=Xmat(:,35);
    


        
    
    

    
 %% bonuses with principal preferences
 
    N=500;  %number of draws over preferences (250 is baseline)   
   % multiplier=[0 1 1.75 2.5 5 10 15 20  40 50 60 65 70 80 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500 525 550 ]';
    %multiplier=[0 10 20 30  50 ]';
    multiplier_fdva=[0 0.0001  0.06 0.08 0.10  0.2 0.3 3.4 3.5 4 4.5 5 5.5 6 7 8.5 10 11.5 13 15 17 19 21 23 25 ]';
    multiplier_fd=[0 0.0001  0.06  0.10 1.9 2 2.5  3 3.5 4 4.5 5 5.5 6 7 8.5 10 11.5 13 15 17 19 21 23 25 ]';

    %multiplier=[ 0 1]';
    mult_fd=multiplier_fd;
    mult_fdva=multiplier_fdva;
    NN=length(multiplier_fd);
    

    stable_school_fd_store=zeros(NN,4);
    stable_teach_fd_store=zeros(NN,4);
    stable_school_fdva_store=zeros(NN,4);
    stable_teach_fdva_store=zeros(NN,4); 
    
    disp('baseline')
    tic
    
 for n=1:NN  
    
     

    stable_school_fd=zeros(N,4);
    stable_teach_fd=zeros(N,4);
    stable_school_fdva=zeros(N,4);
    stable_teach_fdva=zeros(N,4);    
    
    for kk=1:N
[stable_school_fd(kk,:), stable_teach_fd(kk,:), stable_school_fdva(kk,:), stable_teach_fdva(kk,:)] = func_CF_pareto_princ( IJipref,IJipref,IJoutput,IJ_princ, IJVAma,IJclassSize,IJVAm1, IJVAm2,IJVAmean, IJfracDisadv, IJitidx,IJjtidx,IJappyear,thetahat,2015, mult_fd(n,1), mult_fdva(n,1));
    end

    stable_school_fd_store(n,:)=mean( stable_school_fd);
    stable_teach_fd_store(n,:)=mean( stable_teach_fd);
    stable_school_fdva_store(n,:)=mean(stable_school_fdva);
    stable_teach_fdva_store(n,:)=mean(stable_teach_fdva);    
    
    
 end 

 
         savefile= '../../01_data/05_temp/model/cfs_princ';
    save(savefile,'stable_school_fd_store','stable_school_fdva_store','stable_teach_fd_store', 'stable_teach_fdva_store');
toc
   
     
 %% bonuses with VA preferences
 
    N=500;  %number of draws over preferences    
  
   % multiplier=[0 0.0001  0.10 2.5 3 3.5 4 4.5 5 5.5 6 7 8.5 10 11.5 13 15 17 19 21 23 25 ]';
    multiplier=[0 0.0001  0.10 0.85 0.95 2.3 2.5 3 3.5 4 4.5 5 5.5 6 7 8.5 10 11.5 13 15 17 19 21 23 25 ]';

    mult_fd=multiplier;
    mult_fdva=multiplier;
    NN=length(multiplier);
    

    stable_school_fd_store=zeros(NN,4);
    stable_teach_fd_store=zeros(NN,4);
    stable_school_fdva_store=zeros(NN,4);
    stable_teach_fdva_store=zeros(NN,4); 
    
 for n=1:NN  
     n
     
    stable_school_fd=zeros(N,4);
    stable_teach_fd=zeros(N,4);
    stable_school_fdva=zeros(N,4);
    stable_teach_fdva=zeros(N,4);    
    
    for kk=1:N
[stable_school_fd(kk,:), stable_teach_fd(kk,:), stable_school_fdva(kk,:), stable_teach_fdva(kk,:)] = func_CF_pareto_va( IJipref,IJoutput, IJVAma,IJclassSize,IJVAm1, IJVAm2,IJVAmean, IJfracDisadv, IJitidx,IJjtidx,IJappyear,thetahat,2015, mult_fd(n,1), mult_fdva(n,1));
    end

    

    stable_school_fd_store(n,:)=mean( stable_school_fd);
    stable_teach_fd_store(n,:)=mean( stable_teach_fd);
    stable_school_fdva_store(n,:)=mean(stable_school_fdva);
    stable_teach_fdva_store(n,:)=mean(stable_teach_fdva);    
    
    
 end 

 
         savefile= '../../01_data/05_temp/model/cfs_va';
    save(savefile,'stable_school_fd_store','stable_school_fdva_store','stable_teach_fd_store', 'stable_teach_fdva_store');

 